Molecular Dynamics simulations of concentrated aqueous electrolyte solutions 
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Transport properties of concentrated electrolytes have been analyzed using classical molecular dynamics sim- 
ulations with the algorithms and parameters typical of simulations describing complex electrokinetic phenom- 
ena. The electrical conductivity and transport numbers of electrolytes containing monovalent (KCl), divalent 
(MgCl2), a mixture of both (KCl -I- MgCl2), and trivalent (LaCla) cations have been obtained from simulations 
of the electrolytes in electric fields of different magnitude. The results obtained for different simulation parame- 
ters have been discussed and compared with experimental measurements of our own and from the literature. The 
electroosmotic flow of water molecules induced by the ionic current in the different cases has been calculated 
and interpreted with the help of the hydration properties extracted from the simulations. 



I. INTRODUCTION 

Understanding the transport properties of electrolytes in 
aqueous solution is important in a wide range of electrokinetic 
phenomena such as streaming current experiments', ionic 
transport in biological and syntetic nanochannels^-^ or col- 
loidal electrophoresis'*'^. 

Traditionally, electrokinetic and ionic transport phenomena 
have been described using primitive models in which the sol- 
vent is approximated as a continuum of dielectric constant e. 
Although such models provide an accurate description of a 
wide variety of phenomena, they fail for cases in which the 
discrete nature of the solvent plays a fundamental role, or to 
describe phenomena related to specific hydration of ions. 

An alternative theoretical approach to understand elec- 
trokinetic phenomena which circumvents the limitations of 
primitive models is the use of molecular dynamics simula- 
tions of ions in explicit solvent. Due to recent improve- 
ments in algorithms and computer power, all-atom molecular 
dynamics simulatios studies of electrokinetic phenomena in 
some realistic systems are now possible''-^. From such atom- 
istic descriptions, one can elucidate the microscopic mech- 
anisms responsible for macroscopic measurable phenomena. 
For example, Aksimentiev and Schulten studied at atomic 
detail the permeation of individual ions through the trans- 
membrane channel a-Hemolysin with the help of molecu- 
lar dynamics simulations^. This approach has also been em- 
ployed to understand from a microscopic perspective the elec- 
trophoresis of DNA immersed in multivalent ionic solutions 
when an external electric field is applied^, the transport of 
monovalent and divalent ions through polymeric and silica 
nanopores^'"^, and the ionic selectivity of the OmpF porin bi- 
ological nanochannelii. 

The realistic atomic description of such electrokinetic phe- 
nomena involves complex systems containing large number of 
particles. To be able to cope with such systems using molec- 
ular dynamics simulations, the use of algorithms which en- 
hance the computer performance are required. A molecular 
dynamics simulation package which has been proven to be 
very successful in the description of biological molecules and 



elecrokinetic phenomena is NAMD2i2. To accelerate the cal- 
culations of such large systems, in NAMD2 the temperature 
T of the system is controlled by using the Langevin thermo- 
stat instead of the more rigorous but much more computation- 
ally demanding Nose-Hoover thermostat. Also to speed up 
the simulations, NAMD2 is usually employed with a multiple 
time step: a basic timestep for short range interactions and a 
longer one for the evaluation of k-space contribution to the 
long range electrostatic forces in the PME technique. 

In this paper we study, under different conditions, the trans- 
port properties obtained from Molecular Dynamics simula- 
tions of electrolyte solutions using the same algorithms and 
conditions employed in simulations of complex systems. In 
particular, we focus our analysis on the transport proper- 
ties of a monovalent electrolyte (KCl), a divalent electrolyte 
(MgCl2), the mixture of both (KCl + MgCl2) and a trivalent 
electrolyte (LaCla) when an external electric field is applied. 
While the use of the KCl electrolyte is standard in molecular 
dynamics studies and its transport properties have been briefly 
analyzed from this perspective before*", the reliability of such 
simulations methods in dealing with transport properties of 
multivalent electrolytes have not been tested in detail in spite 
of the rich phenomenology that they originat o''^''^'*'* and their 
use in the simulation of electrokinetic phenomena involving 
divalent and trivalent ions^ji^^. 



II. METHODS 
A. Simulation Metliods 

We have studied the transport properties of different elec- 
trolytes by performing molecular dynamics simulations under 
external electric fields of different magnitudes. The systems 
considered were ionic solutions of KCl, MgCl2, a mixture of 
KCl and MgCl2, and LaCls, see Table I. To study the effect of 
the system size on the ionic transport properties of bulk elec- 
trolytes extracted from molecular dynamics simulations, two 
cubic simulation boxes of different size were used, L ~ 4nm 
and L ~ 8nm. 



All simulations have been carried out with the molecular 
dynamics simulation package NAMD2'^, since it is widely 
employed in the simulation of biological macromolecules. 
Water was described using the TIP3P water model as imple- 
mented by the CHARMM force field. The ions were mod- 
eled as charged Lennard-Jones particles with parameters given 
by the CHARMM force field for K+, Mg2+ and for Cr, 
while the Lennard-Jones parameters for La'^^ were taken from 
Ref.'«. 

The initial configuration of the simulation was constructed 
as follows. The ions (K+, Mg^+ and Cl^) were inserted at 
random positions by employing the Autolonize plugin of the 
Visual Molecular Dynamics (VMD) software package'^ in- 
side a cube of preequilibrated TIP3P water molecules, ob- 
tained with the help of the Solvate plugin of VMD (see Fig. [T] 
for a diagram of the system). The faces of the cube are parallel 
to the XY, XZ, and YZ planes. In all simulations we employed 
periodic boundary conditions in all directions. Lennard-Jones 
interactions were computed using a smooth (10-12 A) cut- 
off, as it is customary done in NAMD2 simulations^*^. The 
electrostatic interactions were calculated using the particle- 
mesh Ewald (PME) method with a precision of 10^^ using 
a 128 X 128 x 128 grid and a 12 Acutoff for the real space 
calculation. These are common parameters used to simulate 
complex and large biological macromolecules^. 




FIG. 1: Snapshot of a simulation of the IM MgCb electrolyte with 
simulation box size of L = 3.82 nm. Chloride ions are represented in 
green, magnesium ions in blue, and water molecules in red (oxygen) 
and white (hydrogen). The yellow arrow points in the direction of 
the applied electric field, the green arrow in the direction of the flow 
of chloride ions (anions), and the blue arrow in the direction of the 
flow of magnesium ions (cations). 



For each case, the equilibration procedure consisted of 
50000 steps of energy minimization, a Ins run in the NpT 
ensemble (with p = latm and T ~ 296K) followed by a Ins 
run in the NVT ensemble (T = 296K). Production runs in 
the NVT ensemble were performed in the presence of dif- 
ferent electric fields to induce electromigration of the ions 
in solution. The NpT simulation runs were performed em- 
ploying a combination of the Nose-Hoover constant pressure 



method with the piston fluctuation control implemented using 
Langevin dynamics as implemented in NAMD2 (p ~ latm, 
period of the oscillations of O.lps and relaxation constant of 
0.05 ps). As mentioned above, to speed up calculations the 
NVT runs were carried out by applying Langevin dynamics, 
with parameters (also in the NpT run) T — 296K and a damp- 
ing coefficient of lps~^ (using 0.2ps^^ to test its effect in 
the dynamics for some cases specified later). Langevin forces 
were applied to all atoms except for hydrogens, which ther- 
malize through interactions with the rest of the system. 

In all cases, the equations of motion were solved using a 
multiple time step in order to speed up the simulations. A 
basic time step of 2fs was used for the evaluation of short 
range interactions and a longer time step of 4fs for the calcula- 
tions of the k-space contribution to the long range electrostatic 
forces in the PME technique. 

In the production runs, the instantaneous current induced 
by the external electric field applied along the Z-direction is 
calculated with the help of^ 
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where Zi and qi are the Z-coordinate and the charge of atom 
i, respectively. L is the size of the simulation box and At is 
the time interval employed to record data, which was chosen 
to be At — lOps. The average current is computed by linearly 
fitting the cumulative current that is obtained by integration 
of the instantaneous current. To ensure constistency of the 
results, the current was also computed by counting the flux 
of ions crossing a plane perpendicular to the direction of the 
electric field and located in the center of the simulation box. 
The conductivity k of the solution is defined by; 
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where A — L^ is the cross section area perpendicular to the 
electric field E. In order to obtain the conductivity of the elec- 
trolyte, we have performed simulations at different electric 
fields E and calculated the current / induced by them. In all 
cases, an ohmic behaviour has been observed (i.e. consistent 
with a linear relation between / and E), and a linear fit of the 
data to Eq. (|2]i gives the value of the conductivity k. 

The electromigration of ions in aqueous solutions as a result 
of the application of an external electric field is often accom- 
panied by a net flow of water, the so-called electroosmotic 
flow. The electroosmotic flow has been evaluated by keeping 
track, every At = lOps, of the accumulated number of water 
molecules crossing a plane perpendicular to the electric field. 
Crossings of such plane in the direction of the electric field 
are counted as positive, whereas crossings in the opposite di- 
rection are counted as negative. Hence, the sign of the overall 
flow gives the direction of the electroosmotic flow with re- 
spect to the elecric field direction. As reported before^, the 
simulations might give a drift of the whole system which is 
unphysical since no net force is applied to the system (the 
electrolyte is globally electrically neutral). To avoid such spu- 
rious effects, the computation of the electroosmotic flow has 



TABLE I: Simulation parameters for thie simulations performed of the 
of the Langevin thermostat (rLan), number of K^ ions, number of Mg'^ 
molecules, simulation time, external electric field applied, and potential 
the electric field. 



different electrolytes: cubic box edge length (L), damping constant 
^ ions, number of La'*^ ions, number of CP ions, number of water 
difference between the edges of the cubic box along the direction of 



Electrolyte 


L (nm) 


TLan (pS ^) 


Num. 


Num. 


Num. 


Num. 


Num. 


Simulation 


Electric 


AV (mV) 








K+ 


Mg2+ 


La3+ 


cr 


H2O 


time (ns) 


field (mV/nm) 




IMKCl 


3.88 


1.0 


37 








37 


1907 


21.0 


14.2 


55 




3.88 


1.0 


37 








37 


1907 


9.0 


26.0 


109 




3.88 


1.0 


37 








37 


1907 


9.0 


43.3 


168 




3.88 


1.0 


37 








37 


1907 


9.0 


86.6 


336 


7.82 


1.0 


306 








306 


15774 


9.0 


14.2 


111 




7.82 


1.0 


306 








306 


15774 


9.0 


26.0 


203 




7.82 


1.0 


306 








306 


15774 


9.0 


43.3 


339 




7.82 


1.0 


306 








306 


15774 


9.0 


86.6 


677 


IM MgCb 


3.82 


1.0 





37 





74 


1870 


9.0 


14.2 


54 




3.82 


1.0 





37 





74 


1870 


9.0 


26.0 


99 




3.82 


1.0 





37 





74 


1870 


9.0 


43.3 


165 




3.82 


1.0 





37 





74 


1870 


9.0 


86.6 


331 


7.73 


1.0 





306 





612 


15468 


9.0 


14.2 


110 




7.73 


1.0 





306 





612 


15468 


9.0 


26.0 


201 




7.73 


1.0 





306 





612 


15468 


9.0 


43.3 


335 




7.73 


1.0 





306 





612 


15468 


9.0 


86.6 


669 


7.73 


0.2 





306 





612 


15468 


9.0 


14.2 


110 




7.73 


0.2 





306 





612 


15468 


9.0 


26.0 


201 




7.73 


0.2 





306 





612 


15468 


9.0 


43.3 


335 




7.73 


0.2 





306 





612 


15468 


9.0 


86.6 


669 


O.UMMgCb 


7.83 


1.0 





31 





62 


16293 


9.0 


43.3 


339 


0.33M MgCla 


7.80 


1.0 





93 





186 


16107 


9.0 


43.3 


338 


0.54M MgCla 


7.79 


1.0 





155 





310 


15921 


9.0 


43.3 


337 


IM MgCb 


7.72 


1.0 


306 


306 





918 


14856 


9.0 


14.2 


110 


+ IMKCl 


7.72 


1.0 


306 


306 





918 


14856 


9.0 


26.0 


201 




7.72 


1.0 


306 


306 





918 


14856 


9.0 


43.3 


334 




7.72 


1.0 


306 


306 





918 


14856 


9.0 


86.6 


668 


IM LaCls 


7.71 


1.0 








308 


924 


15154 


13.19 


14.2 


110 




7.71 


1.0 








308 


924 


15154 


10.87 


26.0 


201 




7.71 


1.0 








308 


924 


15154 


10.12 


43.3 


334 




7.71 


1.0 








308 


924 


15154 


10.52 


86.6 


668 



been done in the frame of reference of the center of mass of 
the whole system. 

The electric current flowing in an electrolyte solution is 
caused by the motion of anions and cations moving in oppo- 
site directions under the applied field. The fraction of the total 
current induced by each ion type defines its transport number, 
which is in general a function of the electrolyte concentration. 
The fraction of electrical current carried by cations defines the 
cationic transport number, t_|-, and the fraction of the electrical 
current carried by anions the anionic transport number, t_. A 
completely equivalent quantity to transport numbers is the ra- 
tio of the current induced by anions (anionic current) over the 
current induced by cations (cationic current), which will be 
be refered as the ratio of currents throughout the paper. Due 
to the global electroneutrality of the system, the total electric 



current flowing through an electrolyte as a result of the appli- 
cation of an external electric field is independent of the frame 
of reference in which it is measured. The ratio of currents and 
transport numbers, however, are frame dependent, and their 
computation must be done carefully. A proper account of the 
contribution to the total current from every ion in electrolyte 
solutions is important to describe phenomena in other more 
complex systems in which not only the total electric current 
but also the flow of each type of ion is relevant (e.g. in the 
study of the selectivity of ionic channels). Experimentally, the 
relevant frame of reference in which data of transport num- 
bers and ratio of currents is provided is the frame of reference 
of the moving fluid. Therefore, to facilitate the comparison 
with experimental values, transport numbers and ratio of cur- 
rents will be given in the frame of reference comoving with 



the fluid. 



B. Conductivity measurements 

The electrolyte conductivity measurements were performed 
using a MeterLab CDM210 (R) conductivity meter. Radiome- 
ter Analytical SAS (France). The solutions were prepared us- 
ing water from a Water Purification System Millipore Sim- 
plicity 185. Magnesium Chloride 6-hydrated (MgC12) and 
Potassium Chloride (KCl) from Panreac in all cases following 
ACS specifications. Weighting of the compounds were done 
with a Mettler Toledo AB104-S, in the quantities necessary to 
achieve a IM concentration. The conductivity measurement 
included a stirring of the solution with a magnetic stirrer and 
heater JPSelecta Agimatic-E, with temperature monitoring us- 
ing a laboratory thermometer at 296.0 K. 



III. RESULTS 

A. A test case: Transport properties of IM KCl 

The results obtained from molecular dynamics simulations 
for the electrical properties of the IM KCl electrolyte are 
given in Fig. |2]and Table HI] We have used two different sizes 
for the simulation box, L — 3.88nm and L = 7.82nm, to test 
the dependence of the results on the simulation's box size. As 
can be seen, the values for the ionic conductivity for the two 
different box sizes do not differ significantly, k = 13.4 S/m 
for the smafl box and k — 12.6 S/m for the big one. Both 
results are in good agreement with our measured experimen- 
tal value Kexp — 11.24 ±0.01 S/m, being the value of the 
larger system closer to the experimental result as it should 
be expected. Note that in previous studies of IM KCl bulk 
electrolyte (see Ref.''), it was argued that a 0.2 ps^^ damping 
constant is necessary in order to reproduce correctly transport 
properties. However, our present results were obtained us- 
ing a Langevin thermostat with a damping constant of 1 ps^^, 
which is the typical choice for simulations of complex systems 
in contact with electrolyte (such as protein channels or silica 
nanochanne b^i '°' ^ ' ) . Our results show that with the standard 
simulation parameters employed in complex systems the con- 
ductivity of the electrolyte KCl is correctly reproduced. 

We have also calculated the ratio of the different contribu- 
tions to the total current from anions and cations (equivalent to 
the ratio of transport numbers). As explained in the previous 
section, it is defined as the ratio of the anionic and cationic 
currents in the frame of reference of the moving fluid. The 
results for such ratio of currents are given in Table HI] As 
shown in the table, the values calculated from molecular dy- 
namics simulations are in good agreement with the experi- 
mental value of Id/ Ik = 1-048 (equivalent to a value of the 
cation transport number t+ = 0.488^") for both sizes of the 
cubic simulation box. 
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FIG. 2: Electric current versus the applied electric field for a IM KCl 
electrolyte in simulations using two different simulation box sizes 
L. The points correspond to simulation data computed from Eq.l[T) 
and lines are fits of the data to Eq.l2j (a) L = 3.88 nm (fit gives 
K = 13.4S/m).(b) L = 7.82 nm (fit gives k = 12.6S/m). 



1. Electoosmotic flow 

In addition to ionic currents, we also observe a net flow 
of water molecules. In Fig. JS] the accumulated number of 
water molecules per unit area crossing a plane perpendicular 
to the electric field is represented as a function of simulation 
time. A linear fit of these magnitudes provides the flux of wa- 
ter molecules for every value of the electric field, as shown in 
Table mil As shown in Fig. l3]and Table HIH the electroosmotic 
flow for the IM KCl electrolyte is in the direction opposite to 
the external field, that is, in the du^ection of the flow of chlo- 
ride ions. 



B. Transport properties of IM MgCb 

A similar procedure has been followed to obtain the trans- 
port properties of a IM MgCl2 electrolyte. The electric cur- 
rent induced by the electromigration of ions was obtained for 
different values of the electric field using two different sizes of 



TABLE II: Ionic contributions to the total current, ratio of tlie anionic 
current over tlie cationic current and cationic transport numbers for 
different values of the applied electric field in the case of a IM KCl 
electrolyte. All these quantities are calculated in the frame of refer- 
ence of the moving fluid. 



TABLE III: Flux of water molecules for the different electrolytes and 
different values of the external electric field. 



System Size 


Electric field (mV/nm) 


Ici (nA) 


Ik (nA) 


Ici/Ik 


t+ 




14.2 


1.43 


1.39 


1.03 


0.49 


L = 3.88nm 


25.98 


2.80 


2.59 


1.08 


0.48 




43.3 


4.28 


3.83 


1.12 


0.47 




86.6 


8.42 


8.22 


1.02 


0.49 




14.2 


5.49 


5.14 


1.07 


0.48 


L = 7.82 nm 


25.98 


10.35 


10.11 


1.02 


0.49 




43.3 


17.02 


17.16 


1.00 


0.50 




86.6 


33.57 


33.17 


l.OI 


0.49 
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FIG. 3: Accumulated number of water molecules per unit area cross- 
ing a plane perpendicular to the applied electric field versus simula- 
tion time for the IM KCl electrolyte. Different lines represent the 
flow of water for different values of the electric field. 



a cubic simulation box, L = 3.82 nm and L = 7.73 nm. The 
results of the simulations are given in Figs. |4]and Table HVl 
The values for the ionic conductivity of the two different box 
sizes do not differ significantly, being k — 14.2 S/m for the 
small box and k = 11.9 S/m for the big one. Both results 
agree well with experimental values, being the result from 
the larger simulation box much more accurate. According 
to Phang and Stokes^', ^exp = 11-4 S/m at [MgCl2]=0.9674 
M and T = 298. 15K. A recent critical review^^ provides 
a sUgthly larger value of Kgxp = 11-6 S/m at [MgCl2]=l 
M and T = 25°C. Our own measurements give a value of 
Ke^p = 11.92 S/m. 

The ratio between the different contributions of anions and 
cations to the total current has been evaluated. The results for 
such ratio of currents are given in TablellVI These results em- 
phasize the need for using a large enough simulation box. The 
values for the ratio of currents and transport numbers obtained 
for the cubic simulation box of size L — 3. 82nm exhibit a spu- 
rious dependence on the applied electric field. For the larger 
simulation box, the results are not field dependent, as should 



Electrolyte 


Electric field (mV/nm) 


Water flux (nm ^ns ^) 


IMKCl 


14.2 


0.12 ±0.13 




26.0 


-0.13 ±0.04 




43.3 


-0.43 ±0.14 




86.6 


-0.64 ± 0.07 


IM MgCla 


14.2 


0.76 ± 0.04 




26.0 


1.47 ± 0.05 




43.3 


2.30 ±0.14 




86.6 


5.11 ±0.09 


IM KCl + IM MgCla 


14.2 


0.43 ±0.15 




26.0 


0.88 ± 0.06 




43.3 


1.31 ±0.07 




86.6 


2.82 ± 0.23 


IM LaCl3 


14.2 


-0.25 ±0.11 




26.0 


-0.28 ± 0.25 




43.3 


-0.33 ± 0.08 




86.6 


-0.89 ± 0.24 



TABLE IV: Ionic contributions to the total current, ratio of the an- 
ionic current over the cationic current and cationic transport numbers 
for different values of the applied electric field in the case of a IM 
MgCl2 electrolyte. All these quantities are calculated in the frame of 
reference of the moving fluid. 



System 


Electric 


Ici (nA) 


iMg (nA) 


Ici/lMg 


t+ 


Size 


field (mV/nm) 












14.2 


1.60 


0.72 


2.23 


0.31 


L = 3.82nm 


25.98 


3.09 


2.22 


1.39 


0.42 




43.3 


5.63 


3.56 


1.58 


0.39 




86.6 


10.52 


7.42 


1.42 


0.41 




14.2 


5.69 


4.13 


1.38 


0.42 


L = 7.73 nm 


25.98 


11.31 


8.24 


1.37 


0.42 




43.3 


18.05 


13.66 


1.32 


0.43 




86.6 


35.06 


26.30 


1.33 


0.43 



be expected. Such results show that the anionic contribution to 
the current is significantly larger than the cationic contribution 
(being a factor of 1.4 between them). In our simulations these 
differences in the anionic and cationic contributions to the cur- 
rent can be attributed, to a large extent, to differences in dif- 
fusion coefficients between both ions. In order to disentangle 
the diffusional and correlation contributions to the transport 
number, we have evaluated the translational diffusion coeffi- 
cient of each ion by computing the mean square displacement 
of each ion in a 2ns long NVE simulation run with no exter- 
nal field applied. The results of Dug ~ 0.95 x 10^^ cm^/s 
and Dc\ — 1.69 x 10^^ cm^/s for the diffusion coefficients 
of Mg-^+ and Cl^, respectively, lead to a ratio between the 

diffusion coefficients of _Dci/£'Mg = 1-8. 

Experimentally, both diffusion coefficients and transport 
numbers can be obtained, so we can compare both simulation 
results with experimental data. Using NMR, Struis et al.^^ ob- 
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FIG. 4: Electric current versus applied electric field for a IM MgCb 
electrolyte in simulations using two different simulation box sizes 
L. The points correspond to simulation data computed from Eq.l|T) 
and lines are fits of the data to Eq.l|2j (a) L — 3.82nm size cubic 
simulation box (fit gives k = 14.2S/m). (b) L = 7.73nm (fit gives 
K = 11.9S/m) 



tained Dug = 0.93 x 10"'^ cm^/s and Dci = 2.94 x 10"^ 
cm^/s (at 25 °C and 0.985 mol/Kg concentration), so the ex- 
perimental ratio is Dci/^Mg — 3.16. The simulations re- 
ported here reproduce with good accuracy the diffusion coef- 
ficient for Mg^+ but the diffusion coefficient for Cl~ is greatly 
underestimated. Concerning transport numbers, electrochem- 
ical measurements^ ''^^ give a transport number for cations 
around 0.3, so the experimental ratio between anionic and 
cationic currents is ~ 2.3. The difference between the ratio of 
transport numbers obtained by electrochemical methods and 
the ratio between diffusion coefficients obtained by NMR can 
be interpreted in several ways. First of all, thermodynamic 
arguments^! show that some experimental procedures mixed 
up different reference frames, so the electrochemical results 
have to be interpreted with caution. If the difference between 
both ratios is indeed physical, the difference can be attributed 
to electrokinetic processes (not accounted by diffusion coeffi- 
cients) which are supposed to affect the transport numbers of 
each ion^°. In any case, these electrokinetic processes were 
predicted in the framework of classical, continuum theory of 
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FIG. 5: Diffusion of Mg^^ and CP ions in 1 M aqueous solution as 
obtained from 2ns NVE simulations with a cubic simulation box of 
L = 7.73nm. 



TABLE V: Ionic contributions to the total current, ratio of the an- 
ionic current over the cationic current and cationic transport numbers 
for different values of the applied electric field in the case of a IM 
MgCb electrolyte, L — 7.73nm, and using a damping constant of 
the Langevin thermostat of 0.2 ps~^. All these quantities are calcu- 
lated in the frame of reference of the moving fluid. 




Electric field (V/nm) 


Ici (nA) 


iMg (nA) 


Ici/lMg 


t+ 


14.2 


5.88 


5.25 


1.12 


0.47 


25.98 


11.49 


9.56 


1.20 


0.45 


43.3 


19.82 


15.03 


1.32 


0.43 


86.6 


38.59 


31.15 


1.24 


0.45 



electrolytes and are not clearly observed in our simulations, 
so its molecular origin remains unclear. 

The effect of the damping constant of the Langevin thermo- 
stat was investigated by performing molecular dynamics sim- 
ulations on the same system of size L = 7.73nm and same 
conditions but using a Langevin thermostat with a damping 
constant of TLan = 0.2 ps^^^ (instead of the damping constant 
of 1 ps^^ employed in all other simulations). The electrolyte 
exhibits an ohmic behaviour, with conductivity k = 13.5 S/m, 
which slightly differs from the value obtained before and from 
the experimental value. Nevertheless, the major effect of the 
damping constant is in the different ionic contributions to the 
total current, as can be seen in Fig. [V] In comparison to the 
results of the simulations with TLan = ps~^, it is evident that 
there is a spurious dependence of the ratio of the anionic cur- 
rent over the cationic current on the magnitude of the applied 
electric field. 



1. Electoosmotic flow 

The electroosmotic flow induced by the ionic current was 
also computed. In Fig. |6] the accumulated number of water 
molecules per unit area crossing a plane perpendicular to the 



TABLE VI: Hydration of the different ions 









Number of attached H2O 


Electrolyte 


Ion 


Hydration 


r = lOps 


r = 100 ps 


r = 1000 ps 


IMKCl 


K+ 

cr 


6.4 

7.3 


1.88 
2.41 


0.070 
0.11 


0.0012 
0.0033 


IM MgCl2 


Mg2+ 

cr 


6.0 

7.3 


5.41 
2.99 


5.35 
0.19 


5.31 
0.005 


IM KCI & IM MgCla 


Mg2+ 
K+ 

cr 


6.0 
6.9 

7.0 


5.33 
1.95 
3.04 


5.22 
0.11 
0.25 


5.13 
0.003 
0.009 


IM LaCla 


La-'+ 

cr 


8.4 
6.9 


6.98 
3.50 


4.65 
0.35 


0.27 
0.009 



electric field is represented as a function of simulation time. 
The flux of water molecules for each electric field is given in 
Table |III] It is interesting to note the different direction and 
magnitude of the water flow obtained for IM KCI and for IM 
MgCl2 electrolytes. Indeed, while in the presence of IM KCI 
the direction of the net flow of water is opposite to the di- 
rection of the electric field, in the presence of IM MgCl2 the 
net flow of water goes in the direction of the applied field. 
The magnitude of the net flow of water also differs signifi- 
cantly in both cases, being much larger (almost an order of 
magnitude) for the IM MgCl2 electrolyte. These differences 
can be understood from the hydration properties of the ions 
involved. For each ion, we have computed its hydration, i.e., 
the average number of water molecules in its first coordination 
shell (as defined by the first minimum of the radial distribu- 
tion function, see Figs. [7] and |8]l. We have also computed 
the average number of water molecules which have remained 
a time r in the first coordination shell of each ion to test the 
robustness of the hydration layer, see Table [VT] The values 
given in Table |VI] do not show any dependence on the value 
of the electric field applied in the simulations. The reported 
results for the hydration values of the different ions agree with 
results from previous molecular dynamics studies^^"^^ and Ab 
initio calculations^*^. The results for the average number of 
water molecules that spend a time r in the first coordination 
shell of the different ions are also in agreement with previ- 
ous computational studies, in which residence times of water 
molecules in the first coordination shell of the order of ^ lOps 
were obtained for K+ and Cr and of the order of ^ 10ns 
for Mg^+. Such difference in the robustness of the hydration 
layer of Mg^+ and the other ions is a broadly accepted fact, 
established by ab initio and DFT calculations^^-^'^ as well as 
by NMR, Raman spectroscopy, and X-ray adsorption spec- 
troscopy experiments^^i^Sili. 

Hence, from the analysis of the hydration properties of the 
ions it is possible to understand the difference in the electroos- 
motic flow between the IM KCI and IM MgCl2 electrolytes. 
For the IM KCI electrolyte, the anionic current is slightly 
greater than the cationic current and the hydration of Cr is 
higher than that of K+ so it exhibits a net electroosmotic flow 
in the direction of the flow of Cr . For the IM MgCl2 elec- 
trolyte, although the current of Cr is greater than the current 



of Mg^+, the latter ion is much more effective dragging water 
molecules (its hydration layer is much more robust over time) 
so the net electroosmotic flow is in the direction of the flow of 
Mg-^+ along the applied electric field. 
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FIG. 6: Accumulated number of water molecules per unit area cross- 
ing a plane perpendicular to the applied electric field versus simula- 
tion time for the IM MgCb electrolyte. Different lines represent the 
flow of water for different values of the electric field. 
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FIG. 7; The ion-oxygen and ion-hydrogen radial distribution func- 
tions for the IM KCl electrolyte 
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FIG. 8: The ion-oxygen and ion-hydrogen radial distribution func- 
tions for the IM MgCh electrolyte 



2. Concentration dependence of transport properties ofMgCh 

For the MgCl2 electrolyte we have also tested the depen- 
dency of the conduction properties of the electrolyte on its 
concentration, c. We have performed simulations at differ- 
ent concentrations of MgCl2 for a fixed value of the applied 
electric field E — 0.0433V/nm in a cubic simulation box of 
side L ~ 8nm (see Table |l] for precise values). The summary 
of the results are given in Table IVIII and Figures |9] and [10] 
Here, for each value of the electrolyte concentration the con- 
ductivity has been obtained from a single point (/, E), under 
the assumption that the electrolyte exhibits an ohmic behavior 
(suggested by the studies at IM concentration, see Fig. |4]i. 
As expected, the electrical conductivity, k, of the MgCl2 so- 
lution increases with increasing concentration, as can be seen 



TABLE VII: Total current, electrical conductivity, ionic contribu- 
tions to the total current, ratio of the anionic current over the cationic 
current and cationic transport numbers for different values of the con- 
centration of the MgCl2 electrolyte. 



c(M) 


Itot (nA) 


K(S/m) 


iMg (nA) Ici (nA) lan/Icat t+ 


0.1 


8.0 


3.01 


4.0 


4.0 


1.00 0.5 


0.3 


19.97 


7.55 


9.42 


10.54 


1.12 0.47 


0.5 


28.79 


10.97 


12.90 


15.88 


1.23 0.45 


1.0 


31.70 


12.25 


13.66 


18.05 


1.32 0.43 



in Table IVIll and Fig. |9] In Fig. |9] the conductivity of the 
electrolyte obtained from molecular dynamics simulations for 
each concentration is compared with experimental results — 
and the Kohlrausch's limiting law. The results obtained from 
MD simulations reproduce well the dependency of the elec- 
trolyte conductivity on the concentration, as compared to the 
experimental tendency. Besides, the absolute values of the ex- 
perimental and simulated conductivities agree well with each 
other, especially for high and low salt conditions. In Fig. |9] 
we represent the ratio of the anionic current over the cationic 
current as a function of the MgCl2 concentration. The results 
from MD simulations indicate that the contribution of the an- 
ions becomes more prominent with increasing concentration 
of electrolyte. This tendency agrees with the experimental be- 
haviour for 2 : 1 electrolytes^ii. 
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FIG. 9: Electric conductivity of a MgCl2 electrolyte versus its salt 
concentration for an electric field E = 43.3mV/nm. Blue dots repre- 
sent results obtained from MD simulations, red crosses experimental 
data^', and the green dashed line is Kohlrausch's limiting law^. 



C. Transport properties of mixture composed of IM KCI and 
IM MgCl2 

We have also studied the transport properties of an elec- 
trolyte composed of IM KCl and IM MgCl2. The addition of 
monovalent ions to a system immersed in a multivalent elec- 
trolyte is sometimes used experimentally to screen the elec- 
tric charge of the multivalent ions^ This effect is also used 
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FIG. 10: Ratio of the contribution to the electric current of a MgCl2 
electrolyte by anions and cations versus its bulk concentration when 
an electric field E — 43.3mV/nm is applied. 



in experiment and simulations to determine whether electro- 
static correlations are responsible for a certain macroscopic 
effectli^^^. 

The analysis is similar to the analysis done for the previ- 
ous electrolytes. The electric currents caused by the drift of 
ions have been obtained for four different values of the ap- 
plied electric field. In this case, a single cubic box of size 
L = 7.72 nm was employed. The results for the conductivity 
are summarized in Fig. (TT] which results in a value for the 
conductivity k = 14.8 S/m. Our own experimental measure- 
ments give K ~ 16.83 ± 0.01 S/m. The different contribution 
to the total current of the different ions, as well as the cationic 
transport number are given in Table lVIlU for every value of the 
applied electric field. 
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FIG. 1 1 : MD results for the electric current versus applied electric 
field for an electrolyte composed of IM KCl and IM MgCb in a 
cubic simulation box of L = 7.72nm. The points correspond to 
simulation data computed from Eq.([T]( and lines are fits of the data to 
EqM (fit gives k = 14.8S/m) 



TABLE VIII: Ionic contributions to the total current, ratio of the an- 
ionic current over the cationic current and cationic transport numbers 
for different values of the applied electric field in the case of an elec- 
trolyte composed by IM MgCl2 and IM KCl. All these quantities 
are calculated in the frame of reference of the moving fluid. 



Electric field (V/nm) 


Ici(nA) 


iK(nA) 


iMg(nA) 


i-an/i-cat 


t+ 


14.2 


7.6 


3.27 


2.04 


1.43 


0.41 


25.98 


18.14 


5.12 


4.98 


1.30 


0.43 


43.3 


21.69 


9.29 


7.89 


1.26 


0.44 


86.6 


43.36 


17.11 


15.59 


1.32 


0.43 



1. Electroosmotic flow 

The results for the net flow of water induced by the ionic 
current are given in Fig. [12] in which the accumulated number 
of water molecules per unit area crossing a plane perpendicu- 
lar to the electric field is represented as a function of simula- 
tion time. The flux of water molecules for each electric field 
is given in Table Hill Similarly to the case of the IM MgCl2 
electrolyte, the net electroosmotic flow is in the direction of 
the flow of cations, in spite of being smaller than the net flow 
of anions (see Table [Villi ). We interpret this result along the 
same lines as with the case of IM MgCl2 electrolyte. As it 
is shown in Table IVII the hydration layer of Mg^+ is much 
more robust than the hydration of Cl^ and K+, so the elec- 
troosmotic flow induced by the flow of Mg^+ dominates. 
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FIG. 12: Accumulated number of water molecules per unit area 
crossing a plane perpendicular to the applied electric field versus sim- 
ulation time for an electrolyte composed of IM MgCb and IM KCl. 
Different lines represent the flow of water for different values of the 
electric field. 



D. Transport properties of IM LaCIa 

The lanthanum cation La'^+ is a highly charged and polar- 
izable ion which, for most problems, requires the use of Ab 
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TABLE IX: Ionic contributions to the total current, ratio of the an- 
ionic current over the cationic current and cationic transport numbers 
for different values of the applied electric field in the case of the IM 
LaClg electrolyte, with L = 7.71nm. All these quantities are calcu- 
lated in the frame of reference of the moving fluid. 



Electric field (V/nm) 


Ici (nA) 


iLa (nA) 


Ici/Il. 


t+ 


14.2 


5.51 


5.00 


1.10 


0.48 


25.98 


10.72 


9.89 


1.13 


0.47 


43.3 


17.55 


14.94 


1.17 


0.46 


86.6 


35.00 


31.06 


1.12 


0.47 



initio calculations or polararizable force fields to properly de- 
scribe its interaction with water and ions^'*. However, there 
are systems for which only its electric charge and size are 
relevant to describe certain electrokinetic phenomena^^. For 
such cases, La'^^ can be modeled as a charged Lennard- Jones 
particle and one can obtain the transport properties of a LaCla 
ionic solution by using classical molecular dynamics calcu- 
lations. The Lennard-Jones parameters used here to describe 
the lanthanum cation La'^+ were taken from Ref.^^ 

The electrical transport properties of the IM LaCla elec- 
trolyte in the presence of an external electric field are obtained 
using molecular dynamics simulations which employ a non- 
polarizable force field (described above). The summary of the 
results is given in Fig[T3]and Table llXI The electrolyte shows 
an ohmic behaviour (see FigfTsll. with an electrical conduc- 
tivity of K = 12.8 S/m. Considering the limitations of the 
description, this value is quite satisfactory compared to the 
experimental value of Kexp = 15.3 S/m^^. The ionic contri- 
bution, the ratio of currents and the cationic transport number 
are given in Table |IX] for each value of the external electric 
field. 
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FIG. 13: MD results for the electric current versus applied electric 
field for a IM LaCla electrolyte in a cubic simulation box of L = 
7.71nm. The points correspond to simulation data computed from 
Eq.lUll and lines are fits of the data to Eq.l|2j (fit gives n = 12.8S/m) 



1. Electroosmotic flow 

The results for the net flow of water induced by the ionic 
current are given in Fig. [14] in which the accumulated number 
of water molecules per unit area crossing a plane perpendicu- 
lar to the electric field is represented as a function of simula- 
tion time. The flux of water molecules for each electric field 
is given in Table Hill 

The hydration properties of the IM LaCla electrolyte are 
summarized in Table [Vl] In spite of the limitations of the 
non-polarizable model used to describe La'^^, the values ob- 
tained for the hydration of La^+ are close to the values given 
in the literature (hydration ^^ 8 — 9 and a residence time of 
hydrated water in the first coordination shell of ^ Ins^'*). As 
showed in Fig. [141 the electroosmotic flow obtained for the 
IM LaCla electrolyte is in the direction of the flow of chlo- 
ride ions, opposite to the direction of the applied electric field. 
In this case, the higher hydration of La'^^ is not enough to 
compensate the higher flow of chloride ions versus the flow of 
lanthanum cation La^+. 
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FIG. 14: Accumulated number of water molecules per unit area 
crossing a plane perpendicular to the applied electric field versus 
simulation time for an electrolyte composed of IM LaCla. Different 
lines represent the flow of water for different values of the electric 
field. 



IV. CONCLUSIONS 

Classical molecular dynamics simulations of various elec- 
trolytes (KCl, MgCla, KCl + MgCls, LaClg) in external elec- 
tric fields were performed to study their transport properties. 
We have employed algorithms and parameters typically used 
in simulations of complex electrokinetic phenomena to test 
against experimental data their validity in the description of 
transport properties of electrolytes. The electrical conductiv- 
ity and transport numbers (ionic contributions to the total cur- 
rent) of electrolytes containing monovalent (IM KCl), diva- 
lent (IM MgClz, IM MgCla + IM KCl ), and trivalent (IM 
LaCls) cations were computed from the simulated ion trajec- 
tories. It is shown that in all cases the electrical conductivities 
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FIG. 15: The ion-oxygen and ion-hydrogen radial distribution func- 
tions for the IM LaCla electrolyte 
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electrolyte IM MgCl2, there is a significant discrepancy be- 
tween the transport numbers obtained from MD simulations 
and the experimental values, probably due to the poor result 
for the self difusion coefficient for anions provided by the 
Lennard- Jones parameters used in the simulations. The effect 
of the simulation box size on the electrical transport proper- 
ties was also explored. It is found that a big enough simulation 
box is needed to avoid spurious effects of the simulation. 

The electroosmotic flow of water induced by the ionic flow 
has also been computed from the MD trajectories. It is shown 
that the direction and magnitude of the water flux depends on 
the electrolyte: while the flux is in the direction of the cation 
flow for electrolytes containing Mg^+, it is in the opposite di- 
rection and weaker for all the other cases. These results are in- 
terpreted with the help of the hydration properties of the ions, 
which are calculated and successfully compared with previous 
studies on the subject. 
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